Single cell/nuclei QC and pre-processing

In the chuck below, we will make Seurat objects of all samples (and apply basic filtering based on absolute median deviation to remove outliers):

Description:

• For each object, we keep all genes that are expressed in at least 0 cells, thus we keep all genes as this is required for the removal of ambient RNA (The Cellrange filtered_feature_bc_matrix does not contain genes that have no expression across all cells, but this is often not the case for the raw_feature_bc_matrix, and considering that both matrices will be used for the detection of ambient RNA it is not wise to filter genes with 0 counts in the filtered_feature_bc_matrix. However, this is not an issue if the user chooses to work with the filtered_feature_bc_matrix without including the raw_feature_bc_matrixfrom from the the start).
• Similarly, for now we want to keep all cells that have at least 0 genes expressed. We will later filter out all cells that have not at least x-number of genes expressed after successfully removing ambient RNA. Keep in mind that in reality no cells will be found with 0 genes expressed anyway considering the Cellranger pipeline already removed empty droplets and cells containing <500 UMIs. However, if you choose to work with the raw_feature_bc_matrix, you can actually fish out cell that may truly have <100 genes and <500 UMIs. In that case, we will make sure that we remove droplets (cells) that have <100 genes by default. You can choose to change this parameter, but it is very unlikely to find real cells with <100 genes.

Results:

• From now on you (the user) should keep an eye on the “the number of samples (cells)” in the chuck below. The number of samples in each filtering step will give you an idea how many cell you are losing in each step.

$`ileum_org_pig-G3_IL1_1-filtered_ref_GTF-incl_intr`
An object of class Seurat 
22040 features across 3163 samples within 1 assay 
Active assay: RNA (22040 features, 0 variable features)
 1 layer present: counts

• After filtering the objects based on the absolute median deviation, you can take a close look at the number of samples again and see how many cell were removed when compared to the previous step:
  • Absolute median deviation threshold is set in a way to remove the very extreme cells in case some slipped through the Cellranger filtering.
  • The filtering here will often not remove more than 20-50 cells and will not impact the pre-processing algorithms.
  • This step does have great impact if you choose to work with the raw_feature_bc_matrix. For a good reason too… unless you are in love with empty droplets…
• Again, in the chunk below you should see haw many cells you lost. when compared to the previous step.

$`ileum_org_pig-G3_IL1_1-filtered_ref_GTF-incl_intr`
An object of class Seurat 
22040 features across 3162 samples within 1 assay 
Active assay: RNA (22040 features, 0 variable features)
 1 layer present: counts

QC plots:

Description:

• The number of unique genes detected in each cell:
  • Empty droplets and low quality cells often have very few genes and many 0 counts
  • Cell multiplets (doublets) often have very high gene count.
  • If you work with the filtered counts matrix as a starting point you will not find cells with <500 UMIs.
• User expectation:
  • In general, you should have an expectation of home many genes/UMIs the majority of your cells contain. Most of the time, The UMI counts per cell should be at least ~500 (and ~200 genes), that is considered the low end of what you should expect. However, cell that express extremely low number of genes will have UMI count below ~500. The same is true for a failed experiments with degraded RNAs. Thus the user should have an clear expectation.
• The percentage of reads mapped to the mitochondrial genome:
  • Often, dying cells (or low quality cells) have lots of mitochondrial contamination, possibly due to loss of cytoplasmic RNA from perforated cells. The reasoning is that mitochondria are larger than individual transcript molecules and less likely to escape through tears in the cell membrane.
  • Note that for single-nuclei samples, MT genes are often still present, but at much lower quantities when compared to single-cell data.
  • Mitochondrial QC metrics are calculates with the percentage of counts originating from a set of features.
  • We use the set of all genes starting with “MT-” or “mt-” as mitochondrial genes (which have been inserted/manipulated in ALL MT genes). (The Ensembl database is used to find the MT genes in the species!)

Results:

• Barplots:
  • Representing the distribution and number of gene counts, transcript counts, and percentage of MT reads per cell. Generally, high fraction of MT reads in combination with low gene/transcript counts indicate stressed or dead cells. However, highly prolific stem or cancer-like cell are often rich in MTs, whereas some immune cells such as neutrophiles are in lower end of the spectrum in terms of gene and transcript counts.
  • The dashed red and orange lines represent the theoretical thresholds (Red: 10% MT, 200 nGenes, 550 nUUMIs. Orange: 7.5% MT)

Here we will visualize the distribution of nGenes and nUMIs and the complexity of the cells:

Description:

• Terminology:
  • nFeature_RNA is the number of genes detected in each cell (violin plot).
  • nCount_RNA is the total number of molecules detected within a cell (UMIs) (violin plot).
  • percent.mt is the percentage of mitochondrial UMIs within a cell (violin plot).
  • Cell complexityis the proportion of nGenes raltive to the number of transcripts (RNA-molecules/UMIs).
  • Low complexity cells is a cell that has few genes, but more transcripts.
  • High complexity cells is a cell that more genes and less transcripts.
  • Dominant gene is a gene that has the hiighest expression in a cell relative to the other genes.

Results:

• Violin plots:
  • Low nFeature_RNA for a cell indicates that it may be dead/dying or an empty droplet. It is usually accepted that cells that are poor quality are likely to have low genes and UMIs per cell, and correspond to the data points in the bottom left quadrant of the plot. Good cells will generally exhibit both higher number of genes per cell and higher numbers of UMIs. But be careful here, some cell types are known to have low genes and UMIs.
  • High nCount_RNA and/or nFeature_RNA indicates that the “cell” may in fact be a doublet (or multiplet). But be careful again, some cell types are known to have super high genes and UMIs (i.e. developing cells, cancer, stem-cells).
• Percent-MT vs nGene / Percent-MT vs nUMI plot:
  • Again, for bad cell we expect to see the high percentage MT with low number of genes/UMIs (unless expected otherwise, e.g. cardiac cells).
  • The dashed read line represents a theoretical MT limit of 10% and orange 5%.
• nUMIs vs nGenes plot (also log scaled version):
  • Empty droplets and low quality cells often have very few genes and many “0” genes/UMI counts (Minimum of 500 UMIs if you use the filtered counts matrix).
  • Cell multiplets (doublets) often have very high gene/UMI count.
  • Notice that cell with the highest mitochondrial genes are located in the bottom-left of the plot, indicating that these are truly bad cell considering they contain high MT percentage, but low number of genes/UMIs. Keep in mind that good cells could also have number of MT genes, these cells are often activated cells (i.e. Stem/cancer/cycling/developing cells). Now it could a little troubling if the majority of your high MT cells also contain many genes/UMIs, especially if you do not expect activated cell types in your data. You have to be aware with the type of cell you are working with.
  • In the nUMI vs nGenes plot we can also see the complexity of the dataset with potential errors. Generally, we expect a linear relationship between the number of genes and UMIs in this plot. If cells in the plot are skewed far from the middle line towards the bottom-right this could mean that the experiment has captured very few number of genes and those same genes have been sequenced multiple times. Similarly, if the cells are skewed towards the top-left this could mean that experiment has captured many genes, but they are not sequenced deeply enough (which results in low number of UMIs). Keep in mind that as a result of RNA degradation, a failed experiment could still give us a good plots.
  • The orange dashed lines show the theoretical lower limits of UMIs (550) and genes (200). The yellow dashed lines show the upper and lower quantiles of both number of genes and UMIs (default set between 2-98%). If you choose filter the number of genes and UMIs by user-defined quantiles, then the yellow dashed lines will be set at those numbers.
• Complexity plot (the singular plot below):
  • Complexity can also be visualized by dividing the nGenes by nUMIs for each cell. Higher values indicate that we are getting shallower coverage of more genes, and lower values mean that we are seeing fewer genes overall.
  • Generally, we expect the novelty score to be above 0.80 (indicated with dashed red line).
  • As there is a linear relationship between the number of genes and transcripts for the majority of the cells, it is expected to find the highest probability at a genes/transcript ratio between 0.85-0.95, even for low-complexity cells. A large probability density at anything below 0.8 indicates a failed experiment, as insufficient genes have been captured in proportion to transcript counts.
• Plotting complexity as function of percentage of counts per gene (the first trio plots below):
  • Complexity: Log10 - nUMIs vs nGenes: Again, make another plot as before (log10 - nUMIs vs nGenes plot), but now we fit a linear model and calculate the residual from the middle point for each cell. Now we can clearly visualize which cells have more transcripts vs genes relative to the middle-line (or vice versa). This is indicated by the strength of the blue or red color.
  • Residual vs Density: In this plot, we can see the frequency (density) of cells as a function is the residuals (i.e. deviation from the linear regression). What do you see here? A normal distribution or skewed towards the left or right? This density plot is simply an addition to the linear model and it is expected that you find the largest density around 0, meaning that the majority of the cells were indeed close to the middle line in the linear regression. Just like before, if the graphs is skewed toward the left (negative residuals), this means that you may have an undesired number of low-complexity cells. Likewise, if the graph is skewed towards the right side (positive residuals), then you may many high-complexity cells. Usually, such behavior is related to some sort of failed experiment.
  • Residual vs Percentage Largest genes: In this plot, for each cell, we calculate the gene with the highest percentage of reads and put against the residuals. You can now see a semi-density plot as before, but now with the actual cells. Do you see a pattern here? Do you see any low-complexity cells with a high percentage of a dominant gene? Ideally, we expect a volcano plot centered around the middle. If you have many low-complexity cells with a significantly inflated dominant gene, this might indicate bad cells. But first, we have to find out winch dominant gene they express to make a final conclusion if it is truly related to some biological reason or technical.
• Plotting complexity and examining dominant genes (the final trio plots below):
  • Complexity by Largest Gene: This is the same linear plot as before, but now we visualize the most dominant gene per cell (top 8 out of total in the dataset). Do you see a pattern here or a specific gene that is expressed in low or high complexity cells? Often, the unusual populations in these plots can derive from the activity of a single gene, so now you have a chance figure this out… If these genes happen to be from a mitochondrial origin, then you are likely dealing with dead cells. Likewise, high expression of MALAT1 can also indicate dead cells.
  • Residuals vs Percentage Largest Genes: Same as before, but now with the dominant gene visualized. Which genes are the dominant ones in you off target populations? And does this meet your expectation?
    
  • Residuals vs Percentage Largest Genes + MT %: Same plot, but now with percentage MT. Do the off-target cell populations show high percentage of MT? And do find some correlation between this plot and the dominant genes from before (MT genes)?

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'

Visualize the highy variable genes:

Description:

• We first employ a log normalization method that normalizes the gene expression counts for each cell by the total expression and multiply this by a scale factor of 10,000 and log transform it to stabilize the variance so that larger values do not dominate over smaller ones.

• Then the top 2,500 genes (user-adjustable) with largest standardize variance will be selected automatically and used for the clustering. By selecting highly variable genes, we ensure that house-keeping genes with similar expressions across different cells do not interfere with the unsupervised clustering.

• Note that there is no definitive answer on how many genes one should include for clustering. Thus it is recommended to carefully consider the right number or by trying a range of 2000-5000 variable genes. Generally, low-complexity data (i.e. data with similar cells) will have lower number of variable genes when compared to complex tissue. By default, we will select 2500 highly variable genes.

Results:

• Most Variable Genes:
  • The plot should not looks like a flat line. There should be a clear distribution of the genes with the most highly variable on the top. If there are no highly variable genes in your data, then you are likely dealing with cells that have similar expression levels of the same genes. This can drastically impact the clustering and ambient RNA correction as most downstream steps are computed only on these genes.

Finding variable features for layer counts

Unsupervised clustering: PCA

Description:

• Here, we carry out a linear dimensional reduction using PCA on the scaled data containing the most variable genes (as chosen before):

  • Before putting the highly variable genes into PCA, we will specifically scale those genes so that they have a mean of 0 and a variance of 1. More specifically, For each highly variable gene, a regression model is build using that gene’s expression level across all cells, and then shift the residual to zero and divided it by standard deviation.
  • We are essentially scaling and centering these genes in the dataset, which standardizes the range of expression values. This will regresses out unwanted sources of variation such as technical noise and each gene will contributes similarly to the downstream steps. As a result, the PCA and clustering will not be dominated by the signal of a few highly-expressed genes.

• The reason that we do a PCA is because we simply cannot work with this ultra high dimensional data for the downstream analysis. Instead we can look at which of the most variable genes explains most of the variance in the data.

• Note the PCA is a suitable method, but has some disadvantages:

  • (I still believe that ML would be a better option here, but it is not benchmarked in single-cell data as an alternative to PCA).
  • PCA is very popular, but fails to capture all signals considering that 2D dimensionality reduction cannot capture faint signal.

Results:

• PCA elbow plot: We will find the optimal number of dimensions for downstream analysis, based on:
  • PC where change in percent variation is more than 1%.
  • Or, PCs representing cumulative percent >90% and less than 5% associated with the single PC.
  • With each PC representing combined information of positive and negative correlated gene sets, the top PCs represent strong compression of the total dataset that can be used for further dimensionality reduction and clustering. If you find less than 10 PCs in the elbow plot, your data likely contains too few relevant genes for successful clustering. After all, we are looking for genes that can explain the variance in your cells, and if no such genes can be found (or too few of them) then the downstream steps will have difficulties to find the differences between your cells. However, if you are working with low-complexity data (e.g. similar cell types in different biological states) then you should to find fewer number of PCs.
  • The dashed red line in the plot shows the optimal number PC dimensions.

Final clustering (KNN (SNN)) using the most explained PCs and projection onto UMAP:

Description:

• Original Louvain algorithm:
  • Based on the Euclidean distance in the PCA space, the pipeline carries out a K-nearest neighbor graph (k=20 by default) and refines the edge between any two cells based on shared overlap between their expression patterns using the traditional Jaccard similarity technique.
  • The clusters are then constructed by joining the groups interactively until it fully converges using the Louvain algorithm. Cells are then moved to a non-linear UMAP space for visualization. Here, we find that by using a resolution of 0.8, Louvain does a decent job in identifying sufficient number of communities for both high complexity (porcine ileum, with immune) and low-complexity data (porcine ileum organoid or tissue, no immune).
  • To resolution of 0.8 is set as default, but more research is required to use an adaptive resolution depending on the type of data as this parameter could greatly impact the number preliminary clusters that will be used by the ambient RNA correcting algorithm. As an experimental feature, you are allowed to change the resolution, which will impact the final number of cell clusters (not recommended).
  • In case you end up with too few cluster by the default settings, then it worth to try a resolution >0.8. If you have too few clusters, the ambient RNA correction might become biased (more on this later…).

Results:

• UMAP plots: Number of genes/UMIs/MTs are shown per cell within each cluster:
  • These plots can give us an idea which cells (within clusters) contain low, medium or high number of genes and total UMIs.
  • Note that high number of genes or UMIs is not necessarily a bad thing, but is still highlighted with “purple”
  • In case of a failed experiment (low number of genes) due to significant population of cells with a low RNA content, these plots can show which groups of cells are affected. And if the same clusters contain high percentage MT genes, then it is likely that those cells are low-quality or dead.
  • How many cluster do you observe in total? Do you expect to see more (or less) based on what you already know about your samples?

SoupX Ambient RNA Removal:

Description:

• Terminology:
  • Total contamination is the fraction of ambient RNA (for each gene) estimated from (semi)-empty droplets from the raw counts matrix.
  • Contamination fractions is the the estimated the number transcripts that are contributing to the previously calculated background contamination (from empty droplets) using the expression of several marker genes between cell-clusters where it does not belong. Here it is assumed that the abundance of genes that make up the background are not different between the cells.
• The aim of SoupX is to estimate what fraction of the transcripts are derived as ambient in each droplet and to generate a corrected gene counts matrix with the contaminated counts removed. The algorthm take the following steps:
  • It first calculates the expression profile of the total contamination by estimating the ambient RNA expression from the raw counts-matrix obtained from the CellRanger output. Regardless of the user choice, in this step we always utilize the raw counts as this data contains information of empty droplets, which is an excellent candidate to analyze cell-free droplets that may contain ambient RNA for each gene. (If for some reason the raw counts are not available, SoupX will not run, and the pipeline will continue without the correction of ambient RNA).
  • Once the total contamination has been estimated, SoupX calculates the cell specific contamination level by using the previously found cell clusters. The point here is to find a set of genes in each cell that is known not be expressed in others, and then estimate the contamination fraction using the expression of that given gene in cells it does not belong. If you have too few cluster (due to lack highly variable genes and insufficient number PC dimensions), then it is advised to change number of UMAP dimensions >0.8. Try it out and see if anything improves. Usually, low number of clusters suggests low-complexity data (either biological or technical reasons). Remember that you can also leave out ambient RNA correction when dealing low-complexity if higher UMAP dimensions do not improve anything.
  • With the estimated contamination fraction and the background expression profile at our disposal, the contaminated counts for each gene within each cell can then be corrected.
• Keep in mind that working with the raw counts matrix only will automatically filter cells/droplets that have less than 100 genes. This is a less aggressive filtering when compared to using the filtered counts matrix where Cellranger filtered out cells/droplets that contaied <500 UMIs. For this reason, do not be shocked to find higher ambient RNA contamination when working with the raw counts. The data will be successfully corrected most of the time (unless you are dealing with many degraded RNAs due to failed experiment).

Results:

• Contamination Density plot (First plot):
  

• Representing the probability distribution of estimated contamination fraction. The probability density is calculated by the aggregated set of genes per cells-cluster that have been used to find the contamination fraction. The final contamination fraction is then assigned to the one that occurs with the highest probability (indicated with a vertical red line). A probability distribution with a flat line (no clear peak) indicates failure of SoupX as none of the sets of genes have a preferred contamination fraction, while true estimates are more likely to have the same probability.

• Barplot:

  • Representing the percentage of ambient RNA detected and removed from the data. A contamination level of around 10% is typical for non-nuclei fresh tissue sample.

• Affected Markers Soup Markers (UMAP):

  • Representing the expression of marker candidates that were used for the estimation of the contamination fraction and their level of correction (pre and post Soup correction). By default, top 3 markers are used for plotting, but the user is provided with the option to choose anything between 0-100.

• Affected User Genes (UMAP):

  • Representing the expression of user-defined genes and level of correction. For this to work the user must fill in the genes of interest in the gene.txt file.

• Most Affected Genes (UMAP):

  • Representing the expression of most-affected genes and level of correction. As the most-affected genes are not expressed by many cells, these plots might not be very representative. This option is tuned off by default, but the user may turn it on and select the number plots.

• Colors in the UMAP plots:

  • Red dots - the gene is strongly corrected in the cell
  • White dots - the gene is weakly corrected in the cell
  • Grey dots - Nothing happened to that cell

[1] "Percentage of Ambient RNA (UMI) removed: 17.01 ileum_org_pig-G3_IL1_1-filtered_ref_GTF-incl_intr"

Extensive filtering

Description:

• By default, we automatically filter out cells with a percentage of MT reads above 10%, with gene and transcript counts below the 2% quantile as a lower limit. No filtering is applied on the higher limit as these cells could represent doublets that we want to classify by DoubletFinder instead of removing them by force. However, if the user has prior knowledge on the dataset, filtering can also be applied by user-defined cutoff values for the number of genes, transcripts and percentage of MT reads. Additionally, the user may also adjust the upper and lower quantiles for the parameters mentioned above. It is important to consider that by removing cell in the upper limit by user-defined cutoff values or quantiles, there is a good chance that true doublets are being removed before identified in an algorithmic fashion.

• Note that our default parameters have been performing well in high complexity (porcine ileum and colon tissue, with immune) and low complexity data (porcine ileum organoid) without loosing many cell-types (including some rare immune cells with low number of genes/transcripts).

Results:

• Now take a look in the chunk below… how many cells do you have? Is it significantly lower than before?
  • (The cells are called ‘samples’).
  • (Ignore the inflated number of ‘feautures’. These are the number of genes from two assays in Seurat).

$`ileum_org_pig-G3_IL1_1-filtered_ref_GTF-incl_intr`
An object of class Seurat 
44080 features across 2637 samples within 2 assays 
Active assay: RNA (22040 features, 2500 variable features)
 3 layers present: counts, data, scale.data
 1 other assay present: original.counts
 2 dimensional reductions calculated: pca, umap

Multiplet Removal

Description:

• The DoubletFinder algorithm first generates artificial doublets by randomly merging different cells from the data. In total, one artificial doublet is generated per four cells. We are essentially creating a new dataset where artificial doublets are added to the mix. By doing so, the artificial doublets resemble droplets that may contain two completely different cell types. Arguably, the same process likely occurs as true doublets form in the microfluidic container. DoubletFinder then carries out PCAs using our previously determined optimal number of PC-dimensions (during the unsupervised clustering step). The PC distance matrix is then used to find artificial K-nearest neighbors to define where the true and the artificial doublets share the same neighborhood. This step is recursive and computational expensive and thus we allow the user to choose the number of CPU-cores (4 cores by default). Once complete, we remove the added artificial doublets and the pre-defined proportion of identified real doublets from the dataset.

• Since there is a linear relationship between the percentage of true doublets and number of recovered cells for the “Single Cell 3’ Gene Expression Assay”, we implemented this linearity in the algorithm to ensure the model does not over or underestimate the true doublets in the dataset.

• The relationship between the number of recovered cells and true doublets is provided by 10X (https://kb.10xgenomics.com/hc/en-us/articles/360001378811-What-is-the-maximum-number-of-cells-that-can-be-profiled-). We make a linear regression from this and implement it into DoubletFinder.

• Please increase the number of cores for this step if have many cells. Finding artificial K-nearest neighbors during the parameter sweep can be computationally expensive. The number of cores should NOT exceed the available resources in your system provided to the pipeline. If you do, then DoubletFinder will likely crash!

Results:

• Doublet UMAP plot:
  • Representing the identified doublets within the preliminary cluster. Often, the doublets are located on the edges or between other cell population.
• Violin plots:
  • Representing the scattered distribution of doublets and singlets as a function of gene and transcript counts. Here it is expected that the doublets contain more genes and transcripts when compared to singlets.

ileum_org_pig-G3_IL1_1-filtered_ref_GTF-incl_intr
Doublet Singlet 
     55    2582 

Final filtering and QC plots

Results:

• In the chuck below, we will do another data filtering by taking out the identified muultiplets in the dataset.
  • Now you have can finally compare your results and number of cells before and after filtering. It gives you an idea how aggressively you have filtered your data.
  • Likewise, you can compare the QC plots to what had before. Do you see any improvement here?
  • Remember that you can always change the parameters to your liking an run it again.
• It is advised to try to annotate you cells using the pre-processed data and see if you find your expected cells (in the right proportion). If you think that you are missing some rare cells that contain few genes/transcripts, then try to use the raw counts matrix as a starting point. This will likely expose those cells that have been filtered out by Cellrangers 500 UMI threshold. By using the raw matrix, the pipeline will automatically remove all cells that have <100 genes. By doing so, we ensure that you do not end up many low quality droplets. Even the rarest cells usually do not contain less than 100 genes anyway. Feel free to disagree and change the parameter to your liking. If you do not change anything, then the threshold is always set at 100 genes. NEVER TOUCH THIS PARAMETER IF YOU AEE WORKING WITH THE FILTERED COUNTS MATRIX! You will always lose cells that contain few genes/transcripts.

$`ileum_org_pig-G3_IL1_1-filtered_ref_GTF-incl_intr`
An object of class Seurat 
44080 features across 2582 samples within 2 assays 
Active assay: RNA (22040 features, 2500 variable features)
 3 layers present: counts, data, scale.data
 1 other assay present: original.counts
 2 dimensional reductions calculated: pca, umap

`geom_smooth()` using formula = 'y ~ x'
`geom_smooth()` using formula = 'y ~ x'